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ABSTRACT 

A three node flat shell element with six engineering displacement 
degrees-of-f reedom at each node is developed. The basic 
formulation allows for the arbitrary location of the reference 
surface in which the membrane forces and bending moments are fully 
coupled . 

The well-known, highly accurate, DKT bending element is combined 
with a higher order membrane element in order to obtain a 
consistent formulation. The higher order membrane behavior is 
obtained by the introduction of three additional normal rotational 
degrees-of-f reedom. 

This report presents a summary of the theoretical steps involved 
in the development of the element. The accuracy of the element is 
illustrated by the solution of several standard problems and a 
comparison of results with other thin shell elements. The FORTRAN 
77 listing of the subroutines which form the basic element 
matrices contain less than 300 statements and is presented in 
order to illustrate that the computer implementation of the 
element is relatively simple. 
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INTRODUCTION 

General Background 

The use of composite materials allows for the efficient design of 
many different types of structural systems. One of the major 
advantages of the material is that different stiffness and 
strength properties can be obtained in different directions. 
Therefore, more efficient structures can be obtained since the 
material can be concentrated in the directions of maximum 
stresses . 

Most of the existing finite element programs do not have 
sufficient generality to consider such material properties. Also, 
in the case of thin shell structures very few programs have the 
ability to consider shells in which the bending and membrane 
forces are coupled. In addition, problems associated with the 
modelling of complex shell structures with thin shell elements 
exist since the classical thin shell formulation does not have 
stiffness terms associated with the normal rotational degrees of 
freedom. Therefore, the user of the program is often required to 
add artificial members to a finite element model in order to avoid 
numerical instability in the solution of the finite element 
system. The purpose of this report is to present a new thin shell 
element which is sufficiently robust to solve the above mentioned 
problems . 

Recent Research 

Within the past two years several papers have presented methods 
which introduce a normal rotation in order to improve the membrane 
behavior of plane elements. Carpenter, Stolarski and Belytschko 
present a flat triangular shell element with improved membrane 
interpolation. [1] They introduced normal mid-side displacements 
of the constant strain triangle. The normal displacements are 
eliminated and node rotations are introduced by the use of a cubic 
constraint function along each side of the triangle. A one point 
integration method is used in order to eliminate membrane locking 
within the elements. The element yields very accurate 
displacements; however, the element is rank deficient and is 
unstable for certain geometries. 
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Taylor and Simo applied the same basic approach as presented in 
reference [1] to improve the membrane behavior of quadrilateral 
elements [2] . For many problems excellent displacements and 
stresses are obtained. However, for shell structures such as a 
twisted beam the displacements become very large as the mesh was 
refined. In addition, the flat quadrilateral element does not 
accurately model many common types of shell geometries. Also, the 
DKQ formulation was used to form the bending stiffness which has 
proven to be not as accurate as the DKT formulation. 

Bergan and Felippa have developed a triangular membrane element 
with normal rotational degrees-of-f reedom [3] . The formulation is | 
based on the "free formulation" and uses the continuum-mechanics 
definition of rotation. The element passes the patch test and is i 
stable for all applications. The element produces good values of I 
displacements; however, the values for stresses are poor compared I 
to the values obtained from the Taylor quadrilateral [4] . 

The three elements previously mentioned have not been used for 
thin shells in which the membrane and bending forces are coupled. 
Therefore, one of the purposes of this report is to develop an 
element which has a consistent formulation for both the bending 
and membrane behavior. In addition, the problems with the 
instability associated with normal rotational degree-of-f reedom 
will be studied and a simple technique is suggested in order to 
avoid this problem. 
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BASIC EQUATIONS - ORTHOTROPIC MATERIALS 

The 18 X 18 triangular shell element stiffness matrix for a stiffened 
composite material as shown in figure 1 can be directly calculated 
from the following well-known equation: 



K 



D B dA 



( 1 ) 



The 6x6 D matrix relates the forces to the deformations which are 
associated with a differential element of area dA. Including thermal 
deformations the force-deformation relationship can be expressed by 
the following matrix equation; 



i = D « + f 0 (2) 

where 

= [ mi 1 mx z mt t f 1 1 f X X £i X ] ( 3 ) 

and 

= [ ki 1 kx X ki X Cl 1 Cx X Cl X ] ( 4 ) 

The positive definition of these forces and deformations is 
illustrated in figure 2. 



Normally the matrix D cannot be defined directly for complex 
materials. However, the inverse D~* can normally be easily calculated 
from the basic principles of mechanics or determined experimentally. 
Therefore, the terms for D~‘ are normally specified as input to a 
computer program. The numerical values of D are then evaluated within 
the element stiffness subroutine. Hence, the basic behavior of the 
thin shell, including thermal deformations, is expressed in the 
following form: 
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( 5 ) 
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Figure 1. EXAMPLE OF ANISOTROPIC SHELL 




Figure 2. DEFINITION OF POSITIVE FORCES 



5 



Or in terms of matrix notation: 

£. = D-i f + e_o (5) 

Where dT is the temperature change and ai to 04 are the measured 
thermal expansion coefficients. Hence, the thermal forces indicated 
in equation (2) are calculated from: 

f 0 = - D ip (6) 

Each flexibility term in equation (6) has a direct physical meaning. 
For example, the term Ci x is the curvature kn due to a unit force, 
fix = 1. Whereas, the term Cxi is the strain at the reference 

plane caused by the application of a unit bending moment, mii 1. 
It is apparent that these terms can be best determined experimentally 
for complex composite materials. Also, the values of these terms are 
dependent on the definition of the reference plane which must be 
defined at the same time as the flexibility terms are determined. 
For the special case of constant thickness isotropic shells the 
mid-surface is the logical reference plane and the terms Cij and 
several other flexibility terms are zero in equation (5) . 

In equation (1) the 6x18 B matrix defines the relationship between 
the deformation terms and the node displacements v in the local 



1,2,3 system. 


Or, in 


matrix form: 




A = B V 






(7) 


which can be 


written 


in submatrix form as 




Ip = Bp 


V 




(8a) 


A* ~ B« 


V 




(8b) 



where the "p" and "m” indicate the plate-bending and membrane terms 
respectively. 
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BENDING APPROXIMATION - THE DKT ELEMENT 

The development of the Bp matrix is based on the standard DKT elemen 
[5] . Because the bending and membrane behavior are coupled the DK 
formulation will be summarized here in order that the Bb matrix wil 
be developed with consistent approximations in the next section o 
this report. 

The DKT element is based on the independent expansion of the inplan 
rotations of the reference surface for a 6 node triangle which i 
shown in figure 3. If the local 1 and 2 directions are indicated b 
the local x and y coordinates the finite element approximation i 
written in the following form: 

Px = ai + atx + asy + Ka<x* + asxy + Haey* 

(9) 

pf = a? + as X + a» y + }5ai o x* + ai i xy + J4ai s y* 

The six constants ai to as , ax , can be expressed in terms of the si 
node rotations Px i to Pxs by an inversion of a 6x6 matrix which wil 
produce an equation of the following form: 



ax = H ix (10) 

The same 6x6 matrix, H, will relate the constants a? to ai s , ar , t' 
the node rotations gr . 

From the theory of thin plates the curvature-displacement 
relationships are defined by the following equations: 



kx X 


= Px 


r X 


= at 


+ a<x + aoy 


(llJ 


kf f 


= p» 


# 7 


= at 


+ ai 1 X + a» X y 


(111 


kx f 


M 

GO. 

II 


#7 


+ ?7 #X 






. s 


= as 


+ 


atx + 


asy-t- as -t- atox-i- any 


(IK 



Or, in the following matrix form: 



Cp — Gp a 



( 12 ; 
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(a) BASIC ROTATIONAL UNKNOWNS 

Ljj is the length from I to J 




I M J 



(b) DISPLACEMENTS ALONG TYPICAL SIDE 



Figure 3. DISPLACEMENT APPROXIMATION - DKT ELEMENT 
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Where 



^ — [ ai +33 +33 +34 + 3o + 3i + 



38 +3* + 3i 0 + 3i 1 + 3l « ] 



and 



Gj> 



1 0 X y 0 0 
0 0 0 0 0 0 
0 1 0 X y 1 



0 0 0 0 

1 0 X y 

0 X y 0 



(13 



(14 



The six rotational degrees-of-f reedom which are associated withi 
typical side I-J of the element are shown in figure (3a) . The< 
rotations can be transformed to a n-s reference system which i 
parallel and normal to a typical element side as shown in figure (3t, 
The basic DKT constraints are enforced as follows: 



1. The mid-side rotation 9« ■ is set to the average of the valuj 
at point I and J. Or, 

e.. = ( e.1 + 0BJ ) / 2 (15 

2. The s-displacements , above and below the reference plane, ac 
the normal displacement in the z-direction are forced to be 
cubic functions since the transverse shear strain is in the 
s-direction is set to zero. Therefore, the mid-side normal 
rotation must satisfy the following equation: 

enm = - (Obi + ©bj)/ 4 - 3(wi + wj)/(2Lu) (16 

The constraint specified by equation (15) will force the noritil 
displacement along the element sides, above and below the referenE 
plane, to be linear function. Hence, displacement and slope 
compatibility is satisfied along the sides of all elements. Since : 
attempt is made to set the transverse shear strains to zero within t! 
element the name Discrete Kichoff Triangle was selected as the name ( 
this element. 
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Equation (15) and (16) can be summarized in matrix form as 



es ' 
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’l /2 o' 




'es i' 
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" 1/2 


o' 


Os j 


+ 1 . 5/Li j 
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'Wi ■ 


On 




0 -1/4 




Oh i_ 
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-1/4 


On j 




-1 


1 


Wj 



The relationship between the N-S and X-Y coordinate systems are 



Os 




Cos5 


Sin5 


"ex 


Om 




-Sin5 


Cos5 


ev 


and 










ex 




Cos5 


-Sin5 


es 


ev 




SinS 


Cos5 


Oh 



Equation (17) can now be written in the X-Y system as 



Ox 
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T1 


T2 ' 


ex 1 


+ 


"ti 


T2' 


ex J 


+ 1 . 5/Li J 


’-Ts 


Ts 




’wi" 


Or _ 




T2 


T3 _ 


_ ev 1 




T2 


T3 


ev j_ 




Tc 


-Tc 




Wj_ 



Where 



( 17 ) 



(18a) 



(18b) 



(19) 



T1 = 0.5 Cos* 5 - 0.25 Sin* 5 
T2 * 1.5 Sin5 Cos5 
T3 » 0.5 Sin* 5 - 0.25 Cos* 5 
Ts = 1.5 Lij Sin5 
Tc = 1.5 Lij Cos5 

The six rotational degrees-of-freedom associated with the three 
mid-side nodes of the triangle can be eliminated by the application of 
equation (19) . These transformations can be summarized by a matrix 
equation of the following form: 



e = Tp W 



( 20 ) 



where Tp is a 12x9 matrix and W represents a 9x1 vector of 9x i , 9 y i 
and wi for the three nodes of the triangle. 
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MEMBRANE APPROXIMATION 

The membrane behavior of the triangular shell element is based on th 
basic six node quadratic membrane. If the 12 coefficents are define 
by the 12 x 1 vector b the membrane strains can be expressed in th 
following form: 



= Gm b 



( 21 ) 



As in the case of the DKT plate bending element the three mid-sid 
displacements are rotated to the local N-S coordinate of each side a 
shown in figure (4) . In order to maintain displacement compatibilit 
between element the displacement us is assumed to be linear along eac 
side and the displacement un is a cubic function. These assumption 
can be summarized by the following equations for the displacements a 
the mid-side nodes: 



Us = { Us I + Us j ) / 2 

Un = ( Un I + Un J ) / 2 + Li j { 02 i - 02 j ) / 8 



( 22 ) 



These equations can be written in terms of the X-Y coordinate system 
as 



Ux 




'.5 


o' 




Ux I 




’.5 


o' 




Ux J 




-Sin5 


Sin5 




02 I 




= 




















+ • 125Li J 










Uy 

L 




0 


.5 




Uy I 




0 


.5 




Uy j 




Cos5 


-Cos5 




02 J 



(23) 



The six translational displacements at the midside nodes can now b 
eliminated and three rotational unknowns are introduced at th 
vertices by the direct application of equation (23) . The same basi' 
approach which was used in the DKT formulation is now applied in orde: 
to form the matrix equation of the following form: 



b = T« u (24) 

Therefore, the three membrane strains can be written in terms of th^ 
nine node displacements as 



t_M = Gm Tm u (25) 

It is now possible to evaluate the 24 x 24 element stiffness by the 
direct application of equation (1) . 
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Figure 5 TRIANGULAR 18 DEGREE-OF-FREEDOM SHELL ELEMENT 
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COMPUTER PROGRAM IMPLEMENTATION 



The 18 X 18 triangular element stiffness matrix, for the element showi 
in figure 5, is given by equation (1). Where the the 

strain-displacement matrix B can now be written in terms of th> 
bending and membrane submatrices, or 



B = 



Gp 0 Tp 0 

^ Gm 0. 



G(x,y) T 



( 2 ) 



Since the T matrix is not a function of x and y it is possible t. 
rewrite equation (1) in the following form: 



K 



IT 



I GT D G 



dA 



T 



( 2 '| 



The matrix G is very sparse and contains only the terms 1, x and y 
therefore the integral cab be evaluated directly with a minimum o. 
numerical effort as illustrated by the FORTRAN listing given i: 
Appedix A. In addition, an integral reduction factor can be used o) 
selective terms in order to improve the membrane performance a: 
suggested in reference (3) . 



NUMERICAL EXAMPLES 

In order to illustrate the behavior of the element and to compare th< 
results with other shell elements several examples will be presented 

Cantilever Beam 



The beam shown in figure 6 is idealized by a 1x4 rectangular mesh an< 
is subjected to a load of 40 kips at the tip of the cantilever. Th< 
theoretical displacement at the tip, including shearing deformation 
is 0.3558 inches. 



I 



13 



The Taylor quadrilateral element shell yields a displacement of 0.3467 
inches; or an error of -1.02 percent. Note that the rotations at the 
base of the cantilever are set to zero which is inconsistent with the 
existence of shearing deformations. 

The element presented in this report, TSHELL, was used to model this 
beam with two triangles to form each quadrilateral. The completely 
integrated element produces a displacement of 0.2695 inches, or an 
error of -24.3 percent. With a reduced integration factor of 0.5 the 
tip displacement is 0.3726, or an error of +4.5 percent. 

For all example problems presented in this report a reduced 
integration factor of 0.5 is used. The use of the reduced integration 
factor has the major advantage over one point integration is that a 
"rank deficiency" is not introduced into the element. 



/ 

/ 

/ 

/ 

/ 

/ 

/ 

/ 



E = 30,000 V = 1/4 t * 1.0 




T 
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Figure 6. CANTILEVER BEAM EXAMPLE 
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Twisted Beam 



The twisted beam shown in figure 7 has become a standard test problem 
for thin shell elements [6] . Table 1 summarizes the results obtained 
using four different elements. It is important to note that the 
SHELL element does not converge as the mesh is refined. The reason 
for this unacceptable behavior is because the element is flat and it 
cannot model the twisted surface accurately. The new triangular 
element gives good results for this problem. 



WIDTHS u 
DEPTH = 0.32 
TWIST = 90® 

E =29.0 X 10® 
POISSON RATIO = 0.22 
MESH = I2 x2 
Unit loads at end 

Figure 7. TWISTED BEAM EXAMPLE 



Table 1. Deflection Under Load For Twisted Beam 



Element Type and Mesh 


IN-PLANE 


error 


OUT-OF-PLANE 


error 


Reference Values 


0.005424 


0 


0.001754 


0 


NASTRAM QUAD 4 
NASTRAN QUADS 


0.005386 

0.005413 


-0.7 

-0.2 


0.001727 

0.001750 


-1.5 

-0.2 


SHELL (2x12) 
" (4x24) 


0.005779 

0.006849 


+6.4 

+26.2 


0.001993 

0.002750 


+13.7 

+56.9 


TSHELL (2x12) 
" (4x24) 


0.005390 

0.005399 


-0.6 

-0.5 


0.001717 

0.001735 


-2.1 

-1.1 
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Scordells-Lo Roof 

The shell structure shown in figure 8 is a standard test problem [6] . 
Table 2 summarizes the maximum displacement obtained using different 
elements and meshes. For this problem the SHELL element results are 
very good. However, the TSHELL results appear to converge very 
slowly. However, for even the coarse mesh, the results are of 
acceptable accuracy for normal engineering analysis. 




THICKNESS =0.25 
E =4.32x 10® 

POISSON RATIO* 0.0 
LOADING: 90 /unit ores in z dir. 
u = w = 0 on curved edge 
MESH NxN on quadrant 



Figure 8. SCORDELIS-LO CYLINDRICAL SHELL 



Table 2. Maximum Displacement Of Scordelis-Lo Roof 



Element 


Type and Mesh 


VALUE 


ERROR 


Reference Value 


0.3086 ft. 


0.0 % 


SHELL 


N»4 




+5.2 % 




N=6 




+1.7 % 




N=8 




+0.6 % 




N*10 




+0.2 % 


TSHELL 


N»4 


0.3036 


-1.6 % 




N=8 


0.2982 


-3.4 % 




N=12 


0.2997 


-2.9 % 
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Spherical Shell 

A spherical shell subjected to point loads is shown in figure 9 and is 
also a standard test problem [6] . This structure clearly illustrates 
the weakness of the TSHELL element. With a 12x12 mesh the error ir 
displacement is 17.8 percent. The reason for this very slov 
convergence is that the triangular elements essentially add ril 
reinforcement to the very flexible spherical surface. 

RAOUS » 10.0 




(onquodront) 



Figure 9. SPHERICAL SHELL EXAMPLE 



Table 3. 


Displacement of 


Spherical Shell 






Element 


Type and Mesh 


VALUE 


ERROR 


Reference Value 


0.0925 in. 


0.0 


% 


SHELL 


N«4 




-46.1 


% 




N=8 




-3.3 


% 




N=12 




-0.7 


% 


TSHELL 


N=8 


0.051265 


-44.6 


% 




N=12 


0.075974 


-17.8 


% 
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FINAL REMARKS 

A new anisotropic triangular shell element, with normal rotations at 
the nodes has been developed. The element can be connected directly 
to nodes with beam elements without special consideration. 

The general accuracy of the element has been demonstrated. Care 
should be taken if the element is used to model shells which have 
double curvature. 
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APPENDIX A - FORTRAN LISTING OF TRIANGULAR SHELL ELEMENT 



C 



C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 



C 
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SHELLT 

SUBROUTINE SHELLT (S,T,F,D,C,V,X,Y, AREA , PRES , TEMP ) 
IMPLICIT REAL* 8 (A-H,0-Z) 

INFORMATION CALCULATED 

S = 18 X 18 STIFFNESS MATRIX IN X-Y-Z SYSTEM 

T = 18 X 18 FORCE-TRANSFORMATION MATRIX 

F = 18 X 2 FORCES DUE TO PRESSURE AND TEMP. 

INFORMATION SPECIFIED 

D =6x6 MATERIAL PROPERTY MATRIX 

LOCAL FORCES ARE ASSUMMED TO BE IN THE 
ORDER Mil, M22, M12, Nil, N22 and N12 
C =6x1 MATRIX OF THERMAL EXPANSION TERMS 
V = THE DIRECTION COSINE ARRAY WHICH RELATES THE 
LOCAL 1-2-3 TO THE GLOBAL X-Y-Z SYSTEM 
X, Y, Z GLOBAL NODE COORDINATES 
AREA = ELEMENT AREA 
PRES = AVERAGE SURFACE PRESSURE 
TEMP = AVERAGE TEMPERATURE CHANGE 



WRITTEN BY EDWARD L. WILSON, JAN. - JUNE 1987 



DIMENSION V(4,4),KL(24,3),XY(3),E(6),H(6,6), 

S (18,18) ,T (18,18 ) ,D (6,6) ,X(6) ,Y{6) ,AP (10,12) , 

AM ( 10, 12) , A (20, 18) , AIN (3, 3) , G ( 20 , 20 ) , JP ( 9 ) , 

JM(9) ,TT(20,18) ,F(18,2) ,C(6) 

DATA 

JP /4, 5, 10, 11, 16, 17, 3, 9, 15/, 

JM /I, 2, 7, 8, 13, 14, 6, 12, 18/ 

DATA KL / 

1,1,1, 2,2,2, 3, 3, 3, 3, 3, 3, 4,4,4, 5,5,5, 6, 6, 6, 6, 6, 6, 

1,3,4, 7,9,10, 2, 4, 5, 6, 8, 9, 11,13,14, 17,19,20,12,14,15,16,18, 
1,2,3, 1,2,3, 1,2, 3, 1,2, 3, 1,2,3, 1,2,3, 1,2, 3, 1,2, 3/ 

NT = 24 

CHANGE TO LOCAL COORDINATES SYSTEM 

XC= ( X(l) +X(2) +X(3) ) / 3. 

YC= ( Y(l) +Y(2) +Y(3) ) / 3. 

DO 20 1=1,3 
X(I) = X(I) - XC 

Y(I) = Y(I) - YC 

DO 20 J=l,3 
AIN(I,J) = 0.0 

COMPUTE INTEGRALS 

RED = 0.50*AREA 
AIN (1,1) = AREA 



AIN (2, 2) =-RED*( 


X(l) *X(2) 




X(2)*X(3) 


+ 


X(3)*X(1) 


) 


/ 6. 


AIN (2, 3) = RED*( 


X(l) *Y(1) 


+ 


X(2) *Y(2) 




X(3) *Y(3) 


) 


/ 12 


AIN (3, 3) =-RED*( 


Y(l) *Y(2) 


+ 


Y(2) *Y(3) 




Y(3) *Y(1) 


) 


/ 6. 



AIN(3,2) = AIN(2,3) 



j.y - 



FORTRAN LISTING OF TRIANGULAR SHELL ELEMENT 



C FORM "H" ARRAY FOR GENERAL 6-NODE TRIANGLE 

CALL FORMH(H,A,X,Y) 

C FORM "AP" ARRAY FOR PLATE ELEMENT 

CALL FORMAP(AP,H,X,Y) 

C FORM "AM" ARRAY FOR PLANE ELEMENT 

CALL FORMAM(AM,H,X,Y) 

C FORM 20x20 "G" ARRAY 

DO 160 K=l,20 
DO 160 L=1,K 
160 G(K,L) = 0.0 
C 

DO 200 N=1,NT 

I = KL(N,1) 

K = KL(N,2) 

II = KL(N,3) 

DO 190 M=1,NT 
J = KL(M,1) 

IF (D(I, J) .EQ.0.0) GO TO 190 
L = KL(M,2) 

IF (L.GT.K) GO TO 190 
JJ = KL(M,3) 

G(K,L) = G(K,L) + D(I, J) *AIN(II, JJ) 

190 CONTINUE 
200 CONTINUE 



C 



210 
C 



300 

C 



305 

C 



310 
C 



350 



DO 210 K=l,20 
DO 210 L=1,K 
G(L,K) = G(K,L) 

FORM 20x18 "A" ARRAY 
DO 300 J=l,18 
F(J,1) * 0.0 
F(J,2) = 0.0 
DO 300 1=1,20 
A{I,J) = 0.0 

DO 305 J=l,9 
JJ = JP(J) 

DO 305 1=1,10 
A(I,JJ) = AP(I,J) 



DO 310 J=l,9 
JJ = JM(J) 

DO 310 1=1,10 
A(I+10,JJ) = AM(I,J) 
ROTATE TO GLOBAL X, Y 
DO 350 N=l,6 
NZ = 3*N 
NY = NZ - 1 
NX = NY - 1 
DO 350 1=1,20 
XX = A(I,NX)*V(1,1) + 
YY = A(I,NX)*V(2,1) + 
ZZ = A(I,NX) *V(3,1) + 
A (I, NX) = XX 
A (I, NY) = YY 
A(I,NZ) = ZZ 



Z SYSTEM 



A(I,NY) *V(1,2) 
A(I,NY) *V(2,2) 
A(I,NY) *V(3,2) 



+ A(I,NZ) *V(1,3) 
+ A(I,NZ) *V(2,3) 
+ A(I,NZ) *V(3,3) 






FORTRAN LISTING OF TRIANGULAR SHELL ELEMENT 



C FORM 18x18 ELEMENT STIFFNESS 

DO 380 1=1,20 
DO 380 J=l,18 

CALL DOTP(G(l,I) ,A(1, J) ,SUM,20) 

380 TT(I,J) = SUM 
C 

DO 400 1=1,18 
DO 400 J=1,I 

CALL DOTP(A(l,I) ,TT(1, J) ,SUM,20) 

S(I,J) = SUM 
400 S(J,I) = SUM 

C FORM FORCE-DISPLACEMENT TRANSFORMATION ARRAY 

NO = 0 

DO 500 N=l,3 
XY(1) = 1.0 
XY(2) = X(N) 

XY(3) = Y(N) 

C 

DO 450 K=l,18 
C 

DO 410 1=1,6 
410 E(I) = 0.0 
C 

DO 420 M=1,NT 

I = KL(M,1) 

L = KL(M,2) 

J = KL(M,3) 

420 E(I) = E(I) + XY(J)*A(L,K) 

C 

DO 440 1=1,6 
SUM =0.0 
DO 430 J=l,6 

430 SUM = SUM + D(J,I)*E(J) 

II = NO + I 
440 T(II,K) = SUM 

C 

450 CONTINUE 



C 

500 
C 



600 

610 



NO = NO + 6 

FORM THERMAL FORCES 

IF (TEMP. EQ. 0.0) GO TO 625 
DO 600 1=1,6 

CALL DOTP(D(l,I) ,C(1) ,E(I) ,6) 

CONTINUE 

DO 610 1=1,6 

C(I) = - TEMP*E(I) 

DO 620 1=1,18 
F(I,2) = C(l) *A(1,I) 

F(I,2) = C(2)*A(7,I) 

F(I,2) = C(3)*( A(2,I) + A(6,I) ) 
F(I,2) = C(4)*A(11,I) 

F(I,2) = C(5)*A(17,I) 

F(I,2) = C(6)*( A(12,I) + A(16,I) ) 



620 



FORTRAN LISTING OF TRIANGULAR SHELL ELEMENT 



C CALCULATE PRESSURE FORCES 

625 IF(PRES.EQ.O.O) GO TO 800 
FORCE = - AREA*PRES/3.0 
DO 630 1=1,18,6 
F(I ,1) = FORCE*V(l,3) 

F(I+1,1) = FORCE*V(2,3) 

630 F(I+2,1) = FORCE*V(3,3) 

C 

800 RETURN 

C 

END 

Q FORMH 

SUBROUTINE FORMH (H , B , X , Y) 

IMPLICIT REAL*8 (A-H,0-Z) 

DIMENSION H(6,6),B(6,6),Y(4),X(4) 

C FORM COEFFICIENT MATRIX FOR 6 NODE TRIANGLE 



X(4) = 


( 


X(l) 




X{2) 


) 


/ 


2. 


X(5) = 


{ 


X(2) 




X{3) 


) 


/ 


2. 


X(6) = 


( 


X(3) 


+ 


X(l) 


) 


/ 


2. 


Y(4) = 


( 


YU) 




Y(2) 


) 


/ 


2. 


Y(5) = 


( 


Y(2) 




Y(3) 


) 


/ 


2. 


Y(6) = 


( 


Y(3) 




Y(l) 


) 


/ 


2. 


C 
















DO 100 


1= 


= 1,6 












H(1,I) 


= 


1.0 












H(2,I) 


= 


X(I) 












H(3,I) 


= 


Y(I) 












H(4,I) 


= 


X(I) 


*X(I) / 


2. 






H(5,I) 


= 


X(I) 


*Y(I) 








100 H(6,I) 


S 


Y(I) 


*Y(I) / 


2. 







C INVERT TO FORM H MATRIX - 

DO 200 1=1,6 

DO 200 J=1,I 

SUM = 0.0 
DO 190 K=l,6 

190 SUM = SUM + H(I,K) *H(J,K) 
B(J,I) = SUM 
200 B(I,J) = SUM 
C 

CALL SYMSOL(B,H,6,6,0) 

C 

RETURN 

END 
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FORTRAN LISTING OF TRIANGULAR SHELL ELEMENT 



C FORMAP 

SUBROUTINE FORMAP (AP , H , X, Y) 

IMPLICIT REAL*8 (A-H,0-Z) 

DIMENSION AP { 10 , 12 ) , H { 6 , 6 ) , Y { 4 ) , X { 4 ) , IT { 4 ) 

DATA IT /I, 3, 5,1/ 

C FORM 12 DOF MATRIX 

DO 240 1=2,6 
K = 0 

DO 240 J=l,6 
K = K + 1 
AP(I-1,K) = 0.0 
AP(I+4,K) = H(I,J) 

K = K + 1 

AP(I-1,K) = - H{I,J) 

240 AP(I+4,K) = 0.0 

C ELIMINATION OF 4,5,6 MID-SIDE ROTATIONS 

X(4) = X(l) 

Y(4) = Y(l) 

DO 300 N=l,3 

DX = X(N+1) - X(N) 

DY = Y(N+1) - Y(N) 

XL = DSQRT( DX*DX + DY^DY ) 

S = DY / XL 
C = DX / XL 



T1 


= 


C»C/2. - S*S/4 


T2 


= 


0.75*S»C 


T3 


= 


S*S/2. - C*C/4 


XX 




1.5/XL 


TC 


= 


XX*C 


TS 


= 


XX*S 


11 




IT(N) 


12 


= 


11 + 1 


J1 


= 


IT(N+1) 


J2 


= 


J1 + 1 


L2 


= 


6 + N + N 


LI 


= 


L2 - 1 



I = N + 6 
J = I + 1 
IF (N.EQ.3) J = 7 



DO 300 K=l,10 

AP(K,I1) = AP(K,I1) +AP(K,L1)*T1 
AP(K,I2) = AP(K,I2) + AP(K,L1)*T2 
AP(K,J1) = AP(K,J1) +AP{K,L1)*T1 
AP(K,J2) = AP(K,J2) + AP(K,L1)*T2 
W = - AP(K,L1)*TS + AP(K,L2)*TC 
AP(K,L1) = 0.0 
AP(K,L2) = 0.0 
AP(K,I) = AP(K,I) + W 
AP(K,J) = AP(K,J) - W 
300 CONTINUE 



+ AP(K,L2) *T2 
+ AP(K,L2) *T3 
+ AP{K,L2) »T2 
+ AP{K,L2) *T3 



C 



RETURN 

END 
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FORTRAN LISTING OF TRIANGULAR SHELL ELEMENT 



SUBROUTINE FORMAM { AM , H , X , Y) 

IMPLICIT REAL*8 (A-H,0~Z) 

DIMENSION AM(10,12),H(6,6),Y(4),X(4),IT(4) 
DATA IT /I, 3, 5,1/ 

C FORM 12 DOF MATRIX 

DO 240 1=2,6 
K = 0 

DO 240 J=l,6 
K = K + 1 

AM(I-1,K) = H(I,J) 

AM(I+4,K) = 0.0 
K = K + 1 
AM(I-1,K) = 0.0 

240 AM(I+4,K) = H(I,J) 

C ROTATE NODE 4, 5, 6 MID-SIDE DISPLACEMENTS 

X(4) = X(l) 

Y(4) = Y(l) 

C 

DO 300 N=l,3 

DX = X(N+1) - X(N) 

DY = Y(N+1) - Y(N) 



NY = 2*N + 6 
NX = NY - 1 
IX = IT(N) 
lY = IX + 1 
JX = IT(N+1) 

JY = JX + 1 
I = N + 6 
J = I + 1 
IF (N.EQ.3) J = 7 

C ELIMINATE MID-SIDE DISPLACEMENTS 

DO 260 K=l,10 

TT = .125* ( - AM(K,NX)*DY + AM(K,NY)*DX ) 
AM(K,IX) = AM(K,IX) + AM(K,NX)/2. 

AM(K,IY) = AM(K,IY) + AM(K,NY)/2. 

AM(K,JX) = AM(K,JX) + AM(K,NX)/2. 

AM(K,JY) = AM(K,JY) + AM(K,NY)/2. 

AM(K,NX) = 0.0 
AM(K,NY) = 0.0 
AM(K,I) = AM(K,I) + TT 
AM(K,J) = AM(K,J) - TT 
260 CONTINUE 



FORMAM 



C 

300 CONTINUE 



C 



RETURN 

END 



o o o 
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FORTRAN LISTING OF TRIANGULAR SHELL ELEMENT 

Q LOCALT 

SUBROUTINE LOCALT (XYZ , X , Y , V , AREA , lAXIS ) 

IMPLICIT REAL*8 (A-H,0-Z) 

DIMENSION XYZ(3,3) ,X(6) ,Y(6) ,V(4,4) 

TRANSFORM TO LOCAL COORDINATE SYSTEM 

WRITE (*,3000) XYZ 
3000 FORMAT (3F15.4) 

IF (lAXIS.NE.O) THEN 
DO 100 1=1,3 
100 V(I,4) = 0.0 

I = lABS(IAXIS) 

V(I,4) = lAXIS/I 
END IF 

C DEFINE LOCAL 1,2,3 REFERANCE SYSTEM 

CALL VECTOR (V (1,1) , XYZ (1,1) ,XYZ(1,2) ,XYZ(1,3) , 

XYZ (2,1), XYZ (2,2), XYZ (2,3)) 

CALL VECTOR (V (1,2) , XYZ (1,1) ,XYZ(1,2) ,XYZ(1,3) , 

XYZ(3,1) ,XYZ(3,2) ,XYZ(3,3) ) 

CALL CROSS(V(l,l) ,V(1,2) ,V(1,3) ) 

IF (lAXIS.NE.O) CALL CROSS (V(l,4) ,V(1,3) ,V(1,1) ) 

CALL CROSS (V(l,3) ,V(1,1) ,V(1,2) ) 

C 

DO 5 N=l,3 

X(N) = XYZ(N,1) *V(1,1) + XYZ(N,2) *V(2,1) + XYZ (N , 3 ) * V ( 3 , 1 ) 
5 Y(N) = XYZ(N,1) *V(1,2) + XYZ (N, 2 ) *V ( 2 , 2) + XYZ (N, 3 ) *V ( 3 , 2) 

C CALCULATE AREA OF ELEMENT 

AREA = ( X(2)*Y(3) - X(3)*Y(2) + X(3)*Y(1) 

- X(1)*Y(3) + X(1)*Y(2) - X(2)*Y(1) ) / 2.0 
IF ( AREA. LE. 0.0) THEN 
PAUSE ' ZERO OR NEGATIVE AREA ' 

RETURN 
END IF 

C 

RETURN 

END 
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